Operaciones con datos espaciales

Carga de paquetes

library(tidyverse)
library(sf)
library(leaflet)
library(leaflet.extras)
library(leafem)
library(viridisLite)

Carga de datos

Cantones

cantones <- st_read("cantones_2020_simp_10m.geojson", quiet = TRUE) |> 
  st_transform(4326)

Registros de presencia de felidos

# Carga de datos de félidos de Costa Rica
felidos <-
  st_read(
    "felidos.csv",
    options = c(
      "X_POSSIBLE_NAMES=decimalLongitude", # columna de longitud decimal
      "Y_POSSIBLE_NAMES=decimalLatitude"   # columna de latitud decimal
    ),
    quiet = TRUE
  )

# Asignación del CRS WGS84
st_crs(felidos) <- 4326

Visualizacion de ambos conjuntos de datos

# Factor de color basado en los valores únicos de especies
colores_especies <- colorFactor(
  palette = viridis(length(unique(felidos$species))), 
  domain = felidos$species
)

# Mapa leaflet de cantones y registros de presencia de félidos
leaflet() |>
  setView(
    lng = -84.19452,
    lat = 9.572735,
    zoom = 7
  ) |>  
  addTiles(group = "Mapa general (OpenStreetMap)") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales (ESRI World Imagery)"
  ) |>  
  addPolygons(
    data = cantones,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.5,
    popup = paste(
      paste0("<strong>Código del cantón: </strong>", cantones$cod_canton),
      paste0("<strong>Cantón: </strong>", cantones$canton),
      sep = '<br/>'
    ),
    group = "Cantones"
  ) |>  
  addCircleMarkers(
    data = felidos,
    stroke = F,
    radius = 4,
    fillColor = ~colores_especies(felidos$species),
    fillOpacity = 1.0,
    popup = paste(
      paste0("<strong>Especie: </strong>", felidos$species),
      paste0("<strong>Localidad: </strong>", felidos$locality),
      paste0("<strong>Fecha: </strong>", felidos$eventDate),
      paste0("<strong>Fuente: </strong>", felidos$institutionCode),
      paste0("<a href='", felidos$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),    
    group = "Félidos"
  ) |>
  addScaleBar(
    position = "bottomleft", 
    options = scaleBarOptions(imperial = FALSE)
  ) |>  
  addLegend(
    position = "bottomleft",    
    pal = colores_especies,
    values = felidos$species,
    title = "Especies de félidos",
    group = "Félidos"    
  ) |>  
  addLayersControl(
    baseGroups = c(
      "Mapa general (OpenStreetMap)", 
      "Imágenes satelitales (ESRI World Imagery)"
    ),
    overlayGroups = c("Cantones", "Félidos")
  ) |>
  addResetMapButton() |>
  addSearchOSM() |>
  addMouseCoordinates() |>
  addMiniMap(position = "bottomright") |>
  addFullscreenControl()

Datos vectoriales

Creacion de subconjuntos de espaciales

Registros de presencia de felidos (puntos) en el canton de Sarapiqui

# Polígono del cantón de Sarapiquí
sarapiqui <- filter(cantones, canton == "Sarapiquí")

# Puntos de félidos ubicados dentro del cantón de Sarapiquí
felidos_dentro_sarapiqui <-
  st_filter(
    x = felidos,
    y = sarapiqui, 
    .predicate = st_within
  )

plot(sarapiqui$geometry)
plot(felidos_dentro_sarapiqui, add=TRUE)

Warning in plot.sf(felidos_dentro_sarapiqui, add = TRUE): ignoring all but the
first attribute

Pumas (Puma concolor) en el canton de Osa

# Polígono del cantón de Osa
Osa <- filter(cantones, canton == "Osa")

# Registros de pumas 
pumas <- felidos |> filter(species == "Puma concolor"
                           )
pumas_dentro_osa <-
 st_filter(x = pumas, 
   y = Osa,
   .predicate = st_within
   )
   
  
   
plot(Osa$geometry)
plot(pumas_dentro_osa$geometry, add=TRUE)

Union de datos espaciales

Ejemplo: riqueza de especies de felidos en cantones de Costa Rica

1. Union espacial de felidos y cantones. Esto agregara al conjunto de de datos de felidos una columna correspondiente al codigo de canton en que se ubica cada registro de felidos

felidos_union_cantones <-
  st_join(
     x = felidos,
     y = dplyr::select(cantones, cod_canton),
     join = st_within
  )

2. Conteo de la cantidad de especies de felidos en cada canton (por codigo de canton)

riqueza_especies_felidos_cantones <- 
  felidos_union_cantones |>
  st_drop_geometry() |>
  group_by(cod_canton) |>
  summarise(riqueza_especies_felidos = n_distinct(species, na.rm = TRUE))

3. Unión no espacial de cantones con el dataframe con el conteo de especies en cantones (esto le agrega a cada cantón la cantidad de especies de félidos)

cantones_union_riqueza <-
  left_join(
     x = cantones,
     y = dplyr::select(riqueza_especies_felidos_cantones, cod_canton, riqueza_especies_felidos),
     by = "cod_canton"
  ) |>
  replace_na(list(riqueza_especies_felidos = 0))

4. Generacion del mapa de coropletas

# Paleta de colores de riqueza de especies
colores_riqueza_especies <-
  colorNumeric(
    palette = "Reds",
    domain = cantones_union_riqueza$riqueza_especies_felidos,
    na.color = "transparent"
  )

# Paleta de colores de especies
colores_especies <- colorFactor(
  palette = viridis(length(unique(felidos$species))), 
  domain = felidos$species
)

# Mapa leaflet
leaflet() |>
  setView(
    lng = -84.19452,
    lat = 9.572735,
    zoom = 7) |>
  addTiles(group = "Mapa general (OpenStreetMap)") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales (ESRI World Imagery)"
  ) |> 
  addPolygons(
    data = cantones_union_riqueza,
    fillColor = ~ colores_riqueza_especies(cantones_union_riqueza$riqueza_especies_felidos),
    fillOpacity = 0.8,
    color = "black",
    stroke = TRUE,
    weight = 1.0,
    popup = paste(
      paste("<strong>Cantón:</strong>", cantones_union_riqueza$canton),
      paste("<strong>Riqueza de especies:</strong>", cantones_union_riqueza$riqueza_especies_felidos),
      sep = '<br/>'
    ),
    group = "Riqueza de especies"
  ) |>
  addScaleBar(
    position = "bottomleft", 
    options = scaleBarOptions(imperial = FALSE)
  ) |>    
  addLegend(
    position = "bottomleft",
    pal = colores_riqueza_especies,
    values = cantones_union_riqueza$riqueza_especies_felidos,
    group = "Riqueza de especies",
    title = "Riqueza de especies"
  ) |>
  addCircleMarkers(
    data = felidos,
    stroke = F,
    radius = 4,
    fillColor = ~colores_especies(felidos$species),
    fillOpacity = 1.0,
    popup = paste(
      paste0("<strong>Especie: </strong>", felidos$species),
      paste0("<strong>Localidad: </strong>", felidos$locality),
      paste0("<strong>Fecha: </strong>", felidos$eventDate),
      paste0("<strong>Fuente: </strong>", felidos$institutionCode),
      paste0("<a href='", felidos$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),    
    group = "Registros de presencia"
  ) |>  
  addLegend(
    position = "bottomright",    
    pal = colores_especies,
    values = felidos$species,
    title = "Especies",
    group = "Registros de presencia"    
  ) |>  
  addLayersControl(
    baseGroups = c(
      "Mapa general (OpenStreetMap)", 
      "Imágenes satelitales (ESRI World Imagery)"
    ),
    overlayGroups = c(
      "Riqueza de especies",
      "Registros de presencia"
    )
  ) |>
  addResetMapButton() |>
  addSearchOSM() |>
  addMouseCoordinates() |>
  addFullscreenControl() |>
  hideGroup("Registros de presencia")